#include "StdAfx.h"

#include "Definitions.h"
#include "DDscatMain.h"
#include "DDscatParameters.h"
#include "Complex.h"
#include "Vect3.h"
#include "AbstractTarget.h"
#include "GreenFunctionManager.h"
#include "DielectricManager.h"

void Nearfield(AbstractTarget *currentTarget, DDscatParameters *param, const char *cfle1, const char *cfle2, const char *cfleb1, const char *cfleb2, 
	int myid, real akd, Vect3<real> &ak_tf, int ncomp, real aeff, real wave, Complex *cxadia, DielectricManager *dielec, Complex *cxpol1, Complex *cxpol2, 
	Vect3<Complex> &cxe01_tf, Vect3<Complex> &cxe02_tf, Array4Stacked<Complex> *cxzw)
{
/* **
Given
   NRWORD = length (bytes) of real word
          = 4 for single precision, 8 for double precision
   CMDFFT = FFT method
   CFLB1  = name of output file for nearfield B, for polarization state 1
   CFLB2  = name of output file for nearfield B, for polarization state 2
   CFLE1  = name of output file for nearfield E, for polarization state 1
   CFLE2  = name of output file for nearfield E, for polarization state 2
   MYID
   AKD        = k*d
   AK_TF(1:3)   = (k_x,k_y,k_z)*d in target frame
   IDVOUT     = unit to use for output
   IOCC(1:NAT)=0/1 if site is vacant/occupied
   ICOMP(1:3*NAT) = ICOMP(NAT,3)
              = composition identifier for sites 1-NAT, directions x,y,z
   MXCOMP     = dimensioning info for array CXEPS
   MXPBC      = dimensioning info
   NAT0       = number of occupied sites
   NCOMP      = number of distinct compositions
   NX,NY,NZ   = extent of lattice in x,y,z directions
   IPBC       = 0 to do isolated target
                1 to use periodic boundary conditions
                  in either y or z direction, or both y and z
   IORTH      = 1 if only doing a single polarization
                2 if doing two polarizations for each orientation
   NRFLDB     = 0 to omit calculation of B
                1 to use BSELF to compute B
   GAMMA      = convergence parameter used in PBC calculations
   NAMBIENT   = (real) refractive index of ambient medium
   PYD,PZD    = PBC period/d in y and z directions
   DX(1:3)    = lattice spacing in x,y,z directions/d
                d^3 = DX(1)*DX(2)*DX(3)
   X0(1:3)    = location/d in TF for (I,J,K)=(0,0,0)
   AEFF       = effective radius (phyical units) 
                of target or target unit cell
                aeff = (3*Volume/4*pi)^{1/3}
   WAVE       = wavelength in vacuo (physical units)
   CXADIA(1:3*NAT)=diagonal elements of polarizability tensor
                for dipoles at locations J=1-NAT
   CXEPS(1:MXCOMP)=complex dielectric function for compositions IC=1-MXCOMP
   CXPOL1(1:3*NAT)=polarization (in TF) at locations J=1-NAT
                produced by incident E polarization CXE01_TF
                propagating in direction AK_TF
   CXPOL2(1:3*NAT)=polarization (in TF) at locations J=1-NAT
                produced by incident E polarization CXE02_TF
                propagating in directin AK_TF [used only if IORTH=2]
   CXE01_TF(1:3)= incident polarization vector in TF
   CXE02_TF(1:3)= incident orthogonal polarization vector in TF
   CXZC       = work space needed by ESELF
   CXZW       = work space needed by ESELF

 returns

   CXESCA1(1:3*NAT)=macroscopic electric field at lattice sites 1-NAT
                    generated by the polarization CXPOL1
   CXESCA2(1:3*NAT)=macroscopic electric field at lattice sites 1-NAT
                    generated by the polarization CXPOL2 [only if IORTH=2]
   CXBSCA1(1:3*NAT)=macroscopic=microscopic magnetic field at lattice sites 1-NAT
                    generated by the polarization CXPOL1 [only if NRFLB=1]
   CXBSCA2(1:3*NAT)=macroscopic=microscopic magnetic field at lattice sites 1-NAT
                    generated by the polarization CXPOL2 [only if NRFLB=1 and IORTH=2]
 and writes to file
 if NRFLDB=0
      CFLE1  : P, E_inc, E_sca for incident polarization 1
      CFLE2                    for incident polarization 2 [only if IORTH=2]
 or
 if NRFLDB=1
      CFLEB1 : P, E_inc, E_sca, B_inc, B_sca for incident polarization 1 
      CFLEB2 :                               for incident polarization 2
                                             [only if IORTH=2]

History records of Fortran versions removed.

Copyright (C) 2011,2012 B.T. Draine and P.J. Flatau
Copyright (C) 2011,2012 C++ versions, Choliy V.

This code is covered by the GNU General Public License.
** */
	real dtime = Timeit("Nearfield");

	fprintf(stderr, " >Nearfield x0 = %lf %lf %lf\n", currentTarget->X0().data[0], currentTarget->X0().data[1], currentTarget->X0().data[2]);

	int nx = currentTarget->Nx();
	int ny = currentTarget->Ny();
	int nz = currentTarget->Nz();
	int nxy  = nx * ny;
	int nxyz = nxy * nz;
	int nat3 = 3 * nxyz;

	Vect3<Complex> cxb01_tf, cxb02_tf;
	Complex *cxesca1 = new Complex[nat3];
	Complex *cxbsca1 = NULL;
	if (param->Nrfldb() == NearfieldBMethodBoth)
	{
		cxbsca1 = new Complex[nat3];
	}
	currentTarget->Unreduce(cxadia);

	int j;
	Complex *coff = new Complex[ncomp];
	for(j=0; j<ncomp; ++j)
	{
		coff[j] = Complex((real)3., (real)0.) / (dielec->GetCxeps(j) + (real)2.);
	}
//
	GreenFunctionManager *greenManager = GreenFunctionManager::GetInstance();
	for(int tempIorth=1; tempIorth<=param->Iorth(); ++tempIorth)
	{
		Complex *toUnreduce = (tempIorth == 1) ? cxpol1 : cxpol2;
		currentTarget->Unreduce(toUnreduce);
		greenManager->GetElectric()->Self(param->Cmdfft(), toUnreduce, param->Gamma(), currentTarget->Pyd(), currentTarget->Pzd(), ak_tf, akd, currentTarget->Dx(), cxzw, cxesca1);
//
		if (param->Nrfldb() == NearfieldBMethodBoth)
		{
// allocate CXZG = array that will contain Green function for computing B_sca from P
			bool bAlloc = greenManager->AllocateMagneticBuffer();
			if (bAlloc)
			{
				greenManager->GetMagnetic()->Self(param->Cmdfft(), toUnreduce, param->Gamma(), currentTarget->Pyd(), currentTarget->Pzd(), ak_tf, akd, currentTarget->Dx(), cxzw, cxbsca1);
				if (tempIorth == param->Iorth())
					greenManager->DeallocateMagneticBuffer();
			}
			else
				Wrimsg("Nearfield", "Cannot allocate cxzg");
		}
//
		if (tempIorth == 1)
			cxb01_tf =  cxe02_tf;
		else
			cxb02_tf = -cxe01_tf;
//
// Evaluate incident wave for incident polarization state tempIorth
// for adopted convention E02 = khat x E01
// we have                B01 = khat x E01 = E02
//
// have calculated *microscopic* E field CXESCA1
// now convert to macroscopic E field only changes are within the target material

		short *icompData = currentTarget->Icomp().GetData();			// this loop replaces the commented portion of the code below
		for(int jxyz=0; jxyz<nat3; ++jxyz)
		{
			if (icompData[jxyz] > 0)
				cxesca1[jxyz] *= coff[icompData[jxyz] - 1];
		}
//
//              >>>>> Important Note! <<<<<
// The structure of the WRITE statements below *must* conform to the
// structure of the corresponding READ statements in readnf.f90 
// Any changes must be made in both modules.

		const char *fileName = (param->Nrfldb() == NearfieldBMethodElectric) ? 
			((tempIorth == 1) ? cfle1 : cfle2) :
			((tempIorth == 1) ? cfleb1 : cfleb2);
#ifdef _WIN32
		int modex = O_BINARY|O_WRONLY|O_CREAT;
#else
		int modex = O_WRONLY|O_CREAT;
#endif		
		int iobinFile = open(fileName, modex, S_IWRITE);
		if (iobinFile == -1)
		{
			Errmsg("Nearfield", fileName, "");
			Errmsg("Nearfield", "Error opening iobinFile in Nearfield.", "");
		}
		char Buffer[256];
        sprintf(Buffer, "%s%c", Version(), '\0');
		write(iobinFile, Buffer, strlen(Buffer)+1);
		int ii = VersionNum();
		int nrword = sizeof(real);
		write(iobinFile, &ii, sizeof(int));
		write(iobinFile, &nrword, sizeof(int));
		ii = param->Nrfldb();
		write(iobinFile, &ii, sizeof(int));
		write(iobinFile, &nxyz, sizeof(int));
		write(iobinFile, &currentTarget->Nat0(), sizeof(int));
		write(iobinFile, &nat3, sizeof(int));
		write(iobinFile, &ncomp, sizeof(int));
		write(iobinFile, &nx, sizeof(int));
		write(iobinFile, &ny, sizeof(int));
		write(iobinFile, &nz, sizeof(int));
		currentTarget->X0().Write(iobinFile);
		write(iobinFile, &aeff, sizeof(real));
		write(iobinFile, &param->Nambient(), sizeof(real));
		write(iobinFile, &wave, sizeof(real));
		ak_tf.Write(iobinFile);
		(tempIorth == 1) ? cxe01_tf.Write(iobinFile) : cxe02_tf.Write(iobinFile);
		(tempIorth == 1) ? cxb01_tf.Write(iobinFile) : cxb02_tf.Write(iobinFile);
		write(iobinFile, dielec->GetCxeps(), ncomp * sizeof(Complex));
		for(j=0; j<nxyz; ++j)			// TOREFACT
			currentTarget->Icomp().WriteLine(iobinFile, j);
		write(iobinFile, toUnreduce, nat3 * sizeof(Complex));
		write(iobinFile, cxesca1, nat3 * sizeof(Complex));
		write(iobinFile, cxadia, nat3 * sizeof(Complex));
		if (param->Nrfldb() == NearfieldBMethodBoth)
		{
			write(iobinFile, cxbsca1, nat3 * sizeof(Complex));
		}
		close(iobinFile);
	}
//
	delete [] cxesca1;
	if (param->Nrfldb() == NearfieldBMethodBoth)
		delete [] cxbsca1;
	delete [] coff;

	dtime = Timeit("Nearfield");
}
